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ABSTRACT 

We study the evolution of the interstellar and circumstellar media around massive stars 
(A/ ^ 40-M©) from the main sequence through to the Wolf-Rayet stage by means of radiation- 
hydrodynamic simulations. We use publicly available stellar evolution models to investigate the 
different possible structures that can form in the stellar wind bubbles around Wolf-Rayet stars. 
We find significant differences between models with and without stellar rotation, and between 
models from different authors. More specifically, we find that the main ingredients in the for- 
mation of structures in the Wolf-Rayet wind bubbles are the duration of the Red Supergiant 
(or Luminous Blue Variable) phase, the amount of mass lost, and the wind velocity during this 
phase, in agreement with previous authors. Thermal conduction is also included in our models. 
We find that main-sequence bubbles with thermal conduction are slightly smaller, due to extra 
cooling which reduces the pressure in the hot, shocked bubble, but that thermal conduction does 
not appear to significantly influence the formation of structures in post-main-sequence bubbles. 
Finally, we study the predicted X-ray emission from the models and compare our results with 
observations of the Wolf-Rayet bubbles S308, NGC6888, and RCW58. We find that bubbles 
composed primarily of clumps have reduced X-ray luminosity and very soft spectra, while bubbles 
with shells correspond more closely to observations. 



Subject headings: ISM: bubbles 
— stars: winds, outflows 

Introduction 



Massive stars (M ^ 30M Q ) affect the inter- 
stellar (ISM) and their circumstellar (CSM) media 
due to the combined result of the action of their 
strong stellar winds and ionizing photons. At early 
times, massive stars photoionize the surrounding 
interstellar medium, forming an H II region of 
10 4 K gas around themselves, which then begins to 
expand. The stellar wind from the main-sequence 
(MS) phase interacts with this photoionized re- 
gion, producing a hot, tenuous bubble of shocked 
stellar wind m aterial and a shell of swep t-up pho- 
toionized gas ( Dyson fc Williams] 1997). By the 
end of the main-sequence phase, the massive star 
will be surrounded by a low density, hot bubble of 
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some > 20 pc radius for typical densities in star- 
forming regions, bordered by a thick shell of swept- 
up neutral materia l expanding a t a few km s _1 
(jCappa et al.ll2005l lArthurl l200j . Of course, if 
the underlying ISM was not homogeneous, or if 
the star has a substantial velocity relative to the 
ISM (> 10 km s^ 1 ), the details will vary but the 
general picture remains. 

When the star evolves off the main sequence, 
the mass-loss rate increases substantially as the 
star evolves first towards the red side of the HR 
diagram and then back to the blue. During this 
short-lived phase, which can lead to the star be- 
coming a red or yellow supergiant or a luminous 
blue variable depending on its initial mass, the star 
can deposit half of its mass as a slow, dense, most- 
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likely neutral wind into its immediate environ- 
ment. The most intensive phases of stellar mass 
loss during this period will be episodic and non - 
spherical ejections of material ( Humphreys 2010() . 
In the final Wolf-Rayet (WR) phase of evolution, 
this dense circumstellar medium is photoionized 
by the hot WR star and swept up by the fast WR 
wind. The circumstellar medium of Wolf-Rayet 
stars, therefore, contains the history of the most 
recent, and intense, phases of massive star mass 
loss. 

The slow wind-fast wind interaction has been 
studied numerically by many authors, who all 
show that the interaction leads to instabilities 
(Ray leigh- Taylor and thin shell), which result in 
the corrugation and eventual break up of the 
dense shell of swept-up circumstellar material into 
filaments and clumps. Optical images of the 
nebulae S308, NGC6888 and RCW58 around 
the WR stars WR6, WR136 and WR40, re- 
spectively, show clea r evidence of such a process 
dGruendl et alJl2000l). 

In particular, Garcfa-Segura et all (|l996al lbl) 
were the first to study in detail the purely hy- 
drodynamical evolution of WR bubbles and the 
instabilities which lead to the formation of clumps 
and filaments. They took into account the evolu- 
tion of the stellar wind properties from the main 
sequence through the WR stage for a 35M Q and 
a 6OM star as input to their models. Because 
the winds in the main-sequence and red super- 
giant stages are comparatively steady, these stages 
were studied in one-dimensional spherical sym- 
metry, but the formation of instabilities in the 
slow wind-fast wind interaction at the onset of 
the WR stage was st udied in 2D a xisymmetric 
spherical coordinates. iFrever et al. (|2003l |2006) 
studied the interaction of massive stars (35M Q 
and 60Mq) with their environment with particu- 
lar regard to their effect on the energy balance of 
the interstellar medium. These simulations were 
performed in two-dimensional cylindrical symme- 
try from the beginning of the main sequence stage 
and included the radiative tr ansfer of hydrogen- 
ioniz ing radiation. In this way. lFrever et al.l (|2003 , 
2006) were able to study the evolution of the H II 
region around the stellar wind bubble as well as 
the wind-wind interactions. The stellar wind pa- 
ra meters in this work w e re the same as those used 
bv lGarcia-Segura et al.l (|l996al lbh. 



More recent work has related features in super- 
nova and gamma-ray burst (GRB) afterglow ab- 
sorption spectra to structures in the circumstellar 
medium around the progenitor stars, where the 
progenitor of the GRB event is a WR star and 
the structures are formed during the interaction 
of the WR wind with the dense, slow wind of the 



previous evolutionary stag e (Ramirez- Ruiz et al 
20051: Ivan Marie et alJ li ooa lEldridee et al.l l2006t 



van Marie et al.l 120071) IPwarkadasI (|2007t ) and 



Dwarkadas et al. ( 2010l ) explored the interaction 
of the shock wave from the final supernova ex- 
plosion of the star with the evolved circumstellar 
medium and compared the simulated X-ray emis- 
sion with observations of very young supernova 
remnants. 

The evolving circumstellar media around B 
stars, which do not un d ergo a WR stage, was 
modeled by IChita et all (|2007l l2008t ) where the 
anisotropy of the stellar winds due to fast stellar 
rotation was taken into account, p roducing dis- 
tinctiv e hourglass-shaped structures. Ivan Marie et al 
(|2008j ) also investigated the effect of rapid rota- 
tion on the evolution of a low metallicity 20M Q 
star, the resultant non-spherical density profiles 
in the CSM, and the consequences for the GRB 
afterglow spectra. 

In this paper, we are particularly interested 
in the circumstellar media around Wolf-Rayet 
stars. Many of the > 200 known Galactic Wolf- 
Rayet stars are surrounded by nebulae, visible 
in optical images. Of these, spectros copy iden- 
tifies only 10 as wi nd-blown bubbles 
Esteban et allll992t l T L of which only 2 have been 
detected in X-rays. We expect X-rays to arise 
from the hot, shocked WR wind during the inter- 
action with the ejected slow wind material. The 
nebula S308 arou nd WR6 was detected in soft 
X-ra ys by ROSAT dWriggel ll999h and XMM New- 



ton (jOhu et al.l 120031 ). while NGC 



around 
and 
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Wrigge et al.l 1994, 


1998: Wrigge & Wendkerl2002 


: Wrigee et all 2005; 


Zhekov & Park 2011). Both WR bubbles show 



limb-brightened X-ray morphologies and spectra 
dominated by gas at T ~ 10 6 K, with the X- 
ray emission interior to the optical [O III] emi- 



1 The r emainder are classified as photoionized regions IChi 
Il982f) . 
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sion. This X-ray emission is very soft and will 
be strongly affected by absorption due to neu- 
tral hydrogen along the line of sight. S 308 and 
NGC6888 are two of the closer WR bubbles 
(about 1.5 kpc distant) and the greater distance 
to RCW58 (about 3 kpc) could explain why it has 
not been detected in X-rays. 

The observed X-ray emission from wind-blown 
bubbles has proved difficult to model success- 
fully. There are two main schools of thought: 
on the one hand, if we take a simple hot, 
wind-blown bubbl e , suc h as that described in 
Dvson fc Williams! (| 19971 ). the temperature of the 
hot gas is determined by the stellar wind veloc- 
ity, kT — 3^7tihV^/16, where fi ~ 0.5 is the 
mean particle mass for fully ionized gas. A typical 
Wolf-Rayet wind velocity is ~ 1500 km s -1 , which 
implies temperatures of T > 3 x 10 7 K in the hot 
gas. The densities in the hot bubble are very low 
n < 0.01 cm -3 , and so this model would predict 
low luminosity, hard X-rays. On t h e oth er hand, 



there is the oft-used IWeaver et al.l (|1977l ) model, 
which assumes thermal conduction between the 
hot bubble and the "cold" (10 4 K), swept-up shell. 
Thermal evaporation of dense, cold material into 
the hot bubble raises the density and lowers the 
temperature here. This model predicts high lumi- 
nosity, soft X-ray emission. Both of these models 
are one-dimensional (i.e., radial) in concept and 
do not take into account the break-up of the shell 
due to instabilities nor the shock-clump interac- 
tions which result. The main argument against 
thermal conduction is that even a weak magnetic 
held would restrict the movement of the thermal 
electrons . Diffuse X-ray observations of young star 



Gtid el et aLliauuai). planetary nebu lae ([Unu et al. 
20031 lGuerreroll200a lKastnerll2007l) , as well as the 



clusters (IDunne et al.ll2003t iTownslev et al 
2008ft. 



lanetary nebu lae (|Chu et al 



2003; 



wind-blown bubbles S 308 and NGC 6888 suggest 
that neither model is entirely applicable, since the 
observations indicate that the hot gas has a tem- 
perature 10 6 < T < 10 7 K but the luminosities are 



far low er than those predicted by the I Weaver et al 
(1977) model. 



In this paper, we investigate the possible struc- 
tures that can form around Wolf-Rayet stars by 
means of radiation-hydrodynamical simulations 
taking into account the full time evolution of the 
stellar wind mass-loss rate, terminal velocity and 
stellar ionizing photon rate, as studied by previous 



authors. In addition, we include a time-dependent 
treatment of thermal conduction in our models, 
which has not been done before. We compare re- 
sults obtained from three sets of stellar evolution 



mode ls (jMevnet fc Maederl 120031: lEldridge et al 
2006) for stars with masses between 40M Q and 
60M Q , which include a set of models that take 
into account stellar rotation. By calculating the 
chemical abundances expected in the circumstellar 
nebulae in the post-main-sequence stages of these 
stellar evolution models, we can make compar- 
isons with optical observations. We also calculate 
the expected X-ray emission from our numerical 
simulations under the assumption of collisional 
ionization equilibrium and compare our results 
to the wind-blown bubbles S308, NGC 6888 and 
RCW58. By analyzing the structures formed in 
the circumstellar media of our simulations, the 
chemical abundance history of the CSM and the 
time evolution of the X-ray luminosites, we can 
connect observations to possible mass-loss his- 
tories, and hence discriminate between different 
stellar evolution models. 

In § [5]we describe the numerical method includ- 
ing the treatment of thermal conduction, and in 
§[3] we show the main characteristics of the stellar 
evolution models. Section § [4] presents our one- 
dimensional and two-dimensional results and in 
§ [5] we describe our numerical treatment for the 
simulated luminosity and synthetic spectra, and 
we make comparisons with both optical and X-ray 
observations. Section § [6] discusses our results and 
the effect that differences in the stellar evolution 
models make to the structures and observational 
properties of these objects. Finally, § [7] summa- 
rizes and concludes our findings. 

2. Numerical Method 

In this section we provide a general description 
of the numerical scheme, highlighting new features 
such as thermal conduction. Full details of the 
radiation-hydrodynamics scheme are published in 
the cited references and we describe the main fea- 
tures in Appendix [A] 

2.1. General description of the numerical 
scheme 

Our simulations are carried out in one and two 
dimensions. First, the time-dependent, spherically 
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symmetric system of radiation-hydrodynamic 
equations is solved for the main-sequence and RSG 
phases, during which the stellar winds are reason- 
ably steady and instabilities are relatively unim- 
portant. Then, the ID results are remapped to 
a two-dimensional, cylindrically symmetric (r, z) 
grid, to continue the evolution of the circumstcllar 
medium through the WR stage, where the den- 
sity gradients and wind-wind interactions lead to 
instabilities, which severely disrupt the circum- 
stellar medium. 

The ID spherically symmetric Eulerian hydro- 
dynamic conservation equations are solved us- 
ing a second order, finite volume Godunov-type 
scheme with outflow-only boundary conditions at 
the outer boundary. In 2D, the cylindrically sym- 
metric hydrodynamic conservation equations are 
solved using a hybrid scheme, in which alternate 
time steps are calculated by a Godunov method 
and the van L eer flux-vector splitting method 
van Leerlll982h . This reduces numerical artefacts, 



which can be produced in Godunov-type schemes 
when slow shocks are propagating parallel to one 
of the grid directions. The boundary conditions 
in the 2D scheme are those of no flux through the 
axis of symmetry and outflow-only conditions at 
the outer radial and z boundaries. 

The ID simulations start on a grid of 1000 cells 
with a physical size of 5 pc, which is large enough 
to contain the initial Stromgren sphere. We use 
an expanding grid to avoid large computational 
domains at the beginning of the simulation: when 
the outer shock nears the edge of the grid, more 
(uniform ambient medium) cells are added. Each 
ID calculation has the same physical resolution 
(grid cell size dr = 0.002 pc). By the end of the 
RSG stage, the stellar wind bubble, H II region 
and swept-up neutral shell extend up to 30 pc, 
meaning that the corresponding grids have of or- 
der 6000 cells. The full evolution of the CSM 
around the massive stars is followed in ID from 
the beginning of the main-sequence stage through 
to the WR stage. However, once the star enters 
the WR stage, ID is not adequate to capture the 
full richness of the structures formed during the 
wind-wind interactions and so we remap the re- 
sults onto a 2D grid to continue the calculation in 
cilindrical r — z coordinates. 

The majority of the 2D simulations are carried 
out on a fixed grid of 500 radial by 1000 z-direction 



cells with a corresponding physical size of 10 x 
20 pc 2 (i.e, cell size dr x dz — 0.02 pcx0.02 pc). 
This size is chosen because the wind- wind interac- 
tion takes place within this distance for almost all 
models, and t he radius of observe d WR bubbles 
is R < 10 pc (jGruendl et al.ll2000j ). One model is 
run on a slightly larger grid (750 x 1500 cells) be- 
cause the main wind-wind interaction occurs fur- 
ther from the star, but the spatial resolution is the 
same. 

The transfer of the ionizing radiation is car- 
ried out by the method of short characteristics 
( Raga et all !?)!)!) i adapted for spherical and cylin- 
drical symmetry. We consider only a single point 
source for the ionizing radiation (the star), which 
is positioned on the symmetry axis. The hydrody- 
namics and radiative transfer are coupled through 
the energy equation, where the source term de- 
pends on the photoionization heating and the ra- 
diative cooling rate, and through advection equa- 
tions for the ions and neutrals. The ionization bal- 
ance equation takes into account both photoion- 
ization and collisional ionization together with re- 
combination. This basic code has been exten- 
sively tested and applied to t he evolution of H II 
regions in density g radients ([Hennev et al. 2005; 
Arthur fe Hoarell2006l ). 

The radiative cooling function takes into ac- 
count both cooling in collisionally ionized gas (in 
the hot shocked wind bubble) and in the pho- 
toionized gas in the H II region where collisional 
ionization of hydrogen is not an important con- 
tribution to the cooling rate. Cooling rates are 
generated using the Cloudy photoionization code 



( Ferland et al. 1998h for appropriate stellar pa- 



rameters and used to generate a look-up table. 
The cooling rates extend down to temperatures 
of 100 K but photoionization heating ensures that 
the H II region has a temperature of ~ 10 4 K. 

The stellar wind is inpu t to the grid through the 
Chevalier fe Cleed dl985h thermal wind model, as 



a volume-distributed source of mass and energy. 
This avoids the problem of reverse shocks that can 
occur when there is a decrease in the ram pres- 
sure of a stellar wind input in the form of a region 
of high-density, high- velocity flow near the origin. 
The volume of the stellar wind injection region 
is chosen to be sufficiently small so that the stel- 
lar wind always reaches its terminal velocity and 
shocks outside the source region, but sufficiently 
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large so as to be a reasonable geometrical repre- 
sentation of a sphere on a cylindrical grid (in the 
case of the 2D numerical simulations). 

The variation of the stellar mass-loss rate 
with time is taken fro m publicly available stel- 
lar evolution mode ls ([Mevnet k, Maederl 12003 ; 
Eldridge et al.l 120061 ). The stellar wind velocities 



are calculated using the tabulated stellar parame- 
ters (mass, radius, effective temperature and sur- 
face abundances) and the ionizing photon rates 
are obtained from the appropriate stellar atmo- 
sphere rnodel_s_dis_tnbuted with the Starburst 99 
code ( Leitherer et al. 19991 ). All of this informa- 



tion is stored in tables and interpolated as needed 
using the elapsed simulation time. 

2.2. Thermal Conduction 

We also run a set of models that includes the 
effects of electron thermal conduction to assess the 
importance of this physical process for the dynam- 
ics of stellar wind bubbles and their observable 



prop erties such as X-ray emission ([Weaver et al 



19771 ). Thermal conduction is taken to be a diffu- 
sion process, with the heat flux, q, given by 



-D\7T e , 



(1) 



where D is the diffusion coefficient and T e is 



the el ectron t e mpera ture (jCowie k, McKedll977l ). 
From ISpitzerl (|1962l ) , we have that the electron 
mean free path can be written as 



A e = 2.625 x 10 5 T e 7n e /lnA 



(2) 



where n e is the electron number density and In A 
is the Coulomb Logarithm. For a pure hydrogen 
plasma this can be approximated by 

InA = 9.452 + 3/21nT e - ilnn e (3) 

when T e < 4.2 x 10 5 K and, 



In A = 22.37 + In T e - -lnn e (4) 

when T e > 4.2 x 10 5 K. 

From ISteffen et al.l (|2008l) . the diffusion coeffi- 
cient D is then given by 



D = 7.04 x 10" 



(5) 



where the numerical coeffi cient is the res ult of eval- 
uating physical constants (ISpitzerlll962l) . In spher- 
ical symmetry, the non-linear diffusion equation is 

de dT e 13 / 2 dT e \ 

di = pCv -m = ^d-r{ rD -fr)> (6) 

and we solve this using a Crank-Nicholson-type 



scheme (jPress et al.l 119 92). Thermal conduction 



is also included in the 2D simulations by solving 
the corresponding diffusion equation in cylindrical 
symmetry. The effect of conduction is included as 
an update to the internal energy, e, once the diffu- 
sion equation has been solved for the electron tem- 
perature, and we assume equipartition between 
the electrons and ions. 

When the electron temperature T e is high and 
the electron density n e low, the electron mean free 
path A e becomes large compared to the size scale 
of the hot bubble, and the diffusion approximation 
breaks down. This is th e saturated heat flux limit 
(Co wie fc McKedll977l) that can be described as 



<Zsat 



1.72 x lQ- n T? /2 



(7) 



Saturated conduction is treated by limiting the 
electron mean free path in the following manner 

A = min{0.244Ar, 2.625 x 10 5 T e 2 /n e In A}, (8) 

where Ar is the local grid spacing and the nu- 
merical constant ensures that the conductive 
heat flux can never exceed the saturation limit 



(jSteffen et al.ll2008l ). This formulation means that 
the limiting flux can only be reached at the sharp 
edge of the hot bubble, where |AT e | pa T e . A 
drawback of this method of limiting the electron 
mean free path is that it depends explicitly on the 
numerical resolution Ar and so our results will be 
affected by this. However, since all our ID results 
are run at the same numerical resolution, compar- 
isons between simulations are perfectly valid. 

The temperature found from solving the diffu- 
sion equation for every grid cell is used to update 
the internal energy of that cell once the hydrody- 
namic and radiation steps have been completed. 

3. Stellar Evolution Models 

In this paper we use the publicly available stel- 
lar evolution models for initial stellar masses of 



40 and 60 M Q from iMevnet fc Maederl (j2003l ). 
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hereafter MM2003, and the 40, 45, 50, 55, 
and 60 M@ in i tial mass models from STARS 
dEggletorJll97lllPols et aHll995LlEldridge fe Tout 



2004). Both evolutionary models use the same 
nuclear r eaction rates based on the NACRE 
datab ase (lAngulo et al.l Il999h . and OPAL opac- 
ities (jlglesias fc Roger At low tempera- 



tures, the o pacities come fro m Alexander fc Ferguson! 
(|l994 ) and Ferguson et all (|2005l ). Differences be- 
tween the codes such as assumptions, approxima- 
tions, and solution algorithms, lead to differences 
in the post- main-sequence evolution of the massive 
stars. 

We consider only the Solar metallicity mod- 
els, which for MM2003 means X = 0.705 and 

Y = 0.275, while for STARS this is X = 0.7 and 

Y = 0.28 for Z = 0.02. In addition, MM2003 
provide models including stellar rotation, and so 
we investigate the differences between models with 
and without rotation. The principal effects of an 
initial fast stellar rotation are to increase the main- 
sequence lifetime of the star by about 15-25%, 
and also to extend the lifetimes of the stars in the 
Wolf-Rayet stage (|Mevnet fc Maederll2003h . Fur- 
thermore, the nature of the intense period of post- 
main-sequence mass loss can change, for example, 
the 60M Q MM2003 model without rotation un- 
dergoes a short-lived period of intense mass loss 
that can be interpreted as a Luminous Blue Vari- 
able (LBV) phase, while the corresponding model 
with rotation never becomes an LBV. Models from 
STARS do not include rotation. 

We use the MM2003 models with an initial 
equatorial rotation velocity of 300 km s — 1 . This 
rotation velocity is large enough such that rotation 
effects on the stellar evolution are important, but 
not so large that anisotropics in the stellar wind 
are important in the later stages of stellar evolu- 
tion. MM2003 find that, for an initial rotation 
velocity of 300 km s _1 , during the WR stage the 
ratio of the stellar surface velocity to the break-up 
velocity is very small for the 40 and 60M Q stellar 
models (see their fig. 2) and so we do not expect 
the stellar winds to become strongly anisotropic 
during this stage (jlgnace et al.lll996l ). We do not 
take into account the episodic and non-spherical 
nature of the wind during the m ost intense phases 
of mass loss ([Humphreys! 2010l ) . 

Mass loss is very important in the evolution 
of massive stars. For the mass-loss rates in pre- 




2 3 4 

Time [Myr] 




2 3 
Time [Myr] 

Fig. 1. — Mass-loss rate for top panel 40 Mq 
and bottom panel 60M© models. In each panel, 
solid line — MM2003 no rotation, dashed line - 
MM2003 with rotation, dotted line — STARS. 
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WR evoluti on, MM2003 use t he theoretical mass- 
loss rates of IVink et all (|200ll h whi le STARS uses 
the em pirical mass-loss rates from Ide Jager et al.l 
(|1988l ) . In the WR phase, both ev olution codes use 
the m ass-loss rates proposed by iNugis fc Lamersl 
(l2000h . In Figure [T] we compare the mass-loss 
rates of the MM2003 and STARS models. Al- 
though the details of the main-sequence mass loss 
is very similar in each case, there are appreciable 
differences in the lifetimes of the different post- 
main-sequence stages and the peak mass-loss rate. 
The MM2003 model with rotation has a far longer 
main-sequence lifetime. 

Our radiation-hydrodynamic simulations take 
the mass-loss rates as a function of time directly 
from the stellar evolution models. We also need 
the stellar wind velocities and ionizing photon 
rates for our simulations, and these can be cal- 
culated using the values of the stellar radius, ef- 
fective temperature, and surface abundances from 
the stellar evolution models. The stellar wind ve- 
lo cities for the MM2003 m odels are calculated as 
in Kudritzki et al. ( 19891 ) for the main sequence 
phase, then we use 30 km s _1 for the red/yellow 
supergiant wind (T cff < 8 000 K) or 200 km s" 1 
for the LBV wind (8 000 < T < 25 000 K and 
logM > —3.9), as adopted in the Starburst 99 
code ( Leitherer et al. 19991 and updates). The 
wind velocities in the WR phase (defined by 
surface hydrogen mass fraction X < 0.4 and 
T e g > 25 K) are taken from the tables of 
Smith et al.l (|2002f ). For the STARS models, the 



wind velocities are listed as part of the stellar 
evolut ion tables, using the formulae in lVink et al.l 
(200 1]) for the pre-Wol f-Rayet phase and those of 
Nugis fc Lamersl (|2000l) for the Wolf-Rayet phase. 
We corrected the tabulated values by the appro- 
priate constant factors after noting errors in the 
wind velocities. The resulting wind velocities are 
lower in the RSG phase, even dropping as low as 
5 km s , and are also lower in the WR phase 
compared to the values we adopt for the MM2003 
models. The stellar wind velocities are shown 
in Figure [2] It is apparent that t he WR stel- 
lar wind velocities obtained from the Smith et al 



(|2002l ) tables are extremely high, particularly in 
the final stages of stellar evolution. These veloc- 
ities were t aken from th e ^pp-s pectral type cal- 
ibrations of iPrinia et al. (1990), where spectral 
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Fig. 2. — Stellar wind velocity for top panel 40 M 
and bottom panel 60M© models. In each panel 
the models are, solid line — MM2003 no rotation, 
dashed line - MM2003 with rotation, dotted line 
- STARS. 



type is estimated from the stellar effective tem- 
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perature and surface chemical composition. The 
WR wind veloci t ies ob tained by using the fits of 
Nugis fc Lamersl (j2000h depend on the escape ve- 
locity of the hydrostatic core, the luminosity and 
the surface abundances, and are generally lower. 

The ionizing photon rate was obtained by 
using appropriate stellar atmosphere models 
and integrating over the ionizing flux. For 
both MM2003 and STARS models, we used 
the stellar atmospher e tables from Starburst 99 
( Leitherer et al.lll999l). conc r etely the atmosphere 
models of Pauldrach et al.l < 200lh for the main- 



sequen ce stage, the stellar library of lLeieune et al 



(jl997h 'or the post-main sequenc e stage, and for 
the W R stage we use the models of | rIillier fc Millerl 
(|l998l ). as described by ISmith et all (|2002h . The 
resultant ionizing photon rates as a function of 
time are shown in Figure |3l 

In Figure |4] we plot the evolutionary tracks for 
the 40 M© and 60 M Q MM2003 and STARS mod- 
els together with data points that represent the 
properties of the WR stars WR 6 , WR 40, and 

(|2006l ). These are the 



WR 136 from lHamann et al. 
central stars of the three WR wind-blown bubbles 
S308, RCW58, and NGC6888, respectively. The 
initial masses of these WR stars have been contro- 
versial. We can see that the WN4 star WR 6 falls 
exactly on the evolutionary track of the 40 M Q 
no-rotation model from MM2003 but slightly be- 
low the corresponding model with rotation, while 
it lies between the 40 and 5OM STARS models. 
The WN8 star WR40 seems to fit the 60 M Q mod- 
els from MM2003, but apparently would require a 
much higher mass model from STARS. Finally, the 
WN6 star WR 136 would indicate an initial mass 
lower than 40 M Q in both sets of models. How- 
ever, there are large uncertainties in the luminosi- 
ties of these stars due to variability and distance 
determinations. 

4. Results 

We present results of our one and two-dimensional 
simulations. The ID spherically symmetric sim- 
ulations are used to follow the formation of the 
main-sequence stellar wind bubble and H II re- 
gion and to have a general idea of the evolution 
of the circumstellar medium once the star has 
left the main sequence. Two-dimensional, cylin- 
drically symmetric (r, z) simulations are used to 
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Fig. 3. — Ionizing photon rates for top panel 
40A/q and bottom panel 6OM0 models. In each 
panel the models are, solid line — MM2003 no ro- 
tation, dashed line - MM2003 with rotation, dot- 
ted line — STARS. 
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Fig. 4. — Evolutionary tracks in the T e g- 
luminosity plane. Top panel: MM2003 models. 
Solid lines indicate no stellar rotation, dashed lines 
indicate models with rotation, for 6OM (upper) 
and 40 M (lower) initial masses. Bottom panel: 
STARS models. Solid line - 4OAf , dashed line 
- 5OM , dotted line - 6OM . In each panel, the 
circle, square and triangle represen t the values for 



WR6, WR40 and WR136 from Ha mann et al 



(|2006j ). respectively. The region delimited by the 
thin vertical lines indicates the effective tempera- 
ture range for which the stars would be considered 
yellow supergiants (4800 - 7500 K). 



study the interaction between the fast WR wind 
and the slow, dense RSG wind, where instabili- 
ties are expected to break the swept-up shell into 
clumps. The starting point for the 2D simulations 
are the results of the spherically symmetric sim- 
ulations at the end of the RSG phase, which we 
remap onto the cylindrical grid. 

4.1. ID Results 

4- 1.1. Main-sequence stage 

We start from an initial condition of a cold, uni- 
form, neutral medium with density uq = 100 cm~ 3 
and temperature T = 100 K. The photoionized 
(H II) region forms first. The photoionized gas 
is heated to T ~ 10 4 K and begins to expand 
into the surrounding medium at va > 10 km s . 
The conditions at the ionization front require a 
shock to form in the neutral medium ahead of it, 
which sweeps up and compresses the neutral gas. 
Inside the H II region, the highly supersonic stel- 
lar wind forms a two-shock structure: an outer 
shock sweeps up, heats and compresses the pho- 
toionized gas, and an inner shock brakes and heats 
the stellar wind. The two regions have the same 
pressure but different densities and are separated 
by a contact discontinuity. The outer shocked re- 
gion can cool down quickly to the temperature of 
the H II region due to the high densities in the 
swept-up material, but its temperature is main- 
tained at 10 4 K by photoionization. The shocked 
stellar wind, on the other hand, has a very high 
temperature (T > 10 7 K) and low density and so 
docs not cool efficiently. A hot bubble forms inside 
the H II region, which is in turn surrounded by a 
slowly expanding shell of swept-up neutral mate- 
rial. These structures are well know n from text- 
books (e.g., Dvson fc Williams] 1997 ) and the slow 
time evolution of the stellar parameters does not 
cause any noticeable modifications^. The evolu- 
tion of a stellar wind bubble inside an H II region 
is different to that when ionization is no t taken 
into account (e.g., van Marie et al. 20041 ). Once 
the initial Stromgren region has been formed, the 



2 Two-dimensional simulations of the early stages of stel- 
lar wind bubbles inside H II regions show more compli- 
cated structures. Rapid cooling in the material swept 
up by the outer wind shock forms clumps, which then 
trigger the shadowing instability, leading to complex low- 
velocity (< 20 km s —1 ) kinematics in the photo ionized gas 
jFrever et alJl2003L 120061 : lArthur fc Hoardl2006f ) . 
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pressure in the H II region is much higher than 
that in the surrounding ISM. Moreover, although 
to begin with the pressure in the hot, shocked, 
stellar wind bubble is higher than that in the pho- 
toionized gas, pressure equilibrium with the H II 
is reached after a short time and, thereafter, the 
pressure across the bubble structure is regulated 
by both the H II region and the stellar wind to- 
gether. 

To illustrate the main-sequence evolution, in 
Figure[5]we show the results from a 4OM STARS 
model for simulations with and without thermal 
conduction after t = 10 6 yrs. The hot, shocked 
wind bubble, H II region, and neutral shell are 
cleary identified. The bubble radius in the simu- 
lation with thermal conduction is slightly smaller 
than that without conduction due to increased ra- 
diative losses in the hot bubble. These enhanced 
radiative losses in the hot, stellar wind bubble re- 
duce the pressure here and enable the H II region 
to expand back towards the star, and the phoion- 
ized gas has slightly lower density in the thermal 
conduction simulations that in those without con- 
duction as a result. 

For our particular models, we find that by the 
end of MS stage, the radius of the region af- 
fected by the massive star is 25-30 pc (see Ta- 
bic [I]) and the mass of swept-up neutral material 
is > 10 5 M Q . The exact size of the bubble depends 
mainly on the length of the main-sequence stage 
(e.g., models with rotation have longer main se- 
quences than those without) and the stellar wind 
velocity. Higher stellar wind velocities lead to 
higher pressures in the shocked wind bubbles and 
therefore faster expansion speeds (e.g., the STARS 
models have higher main-sequence stellar wind ve- 
locities, see Fig. [2]). We remark that bubbles 
expanding into denser media will, of course, be 
smaller, since on dimensional grounds the exter- 
nal radius of a stellar wind bubble in an uniform 
medium of density n a is R ~ (M M /V2/n ) 1 / 5 t 3 / 5 
where My/ is the mass-l oss rate and Vw the stella r 
wind velocity (see e.g. Dyson fc Williamsl [1997) . 
while that of an expanding H II region follows 
Ri ~ (S,/a B n 2 ) 3 / 7 4{ 7 t 4 / 7 , where «b and en are 
the recombination coefficient and the sound speed 
in the ionized gas, respectively, and S 1 * is the stel- 
lar ionizing photon rate. Both of these rela tions 
depend inversely on the density (jSpitzer 119781 ). 



4-. 1.2. Post-main-sequence stage 

During the RSG phase the star loses a large 
amount of its mass in a slow, dense wind (see 
Fig. [1} . This wind is subsonic with respect to the 
hot wind bubble left by the main-sequence stellar 
wind and so piles up at the bottom of a r -2 density 
gradient into a thin dense shell without forming a 
two-shock pattern. 

In Figure [5] we show the density, temperature 
and internal energy (i.e. thermal pressure, since 
e = p/[l ~ 1]) as a function of radius at the end of 
RSG phase for the 40 M Q models. We compare the 
results using the MM2003 stellar evolution model 
(without rotation) with those of STARS. For the 
simulations without thermal conduction, we see an 
inner region of RSG wind, which has a r~ 2 den- 
sity profile and temperature T < 10 K, with a 
high central density and bounded by a dense shell 
of partially ionized material. For the MM2003 
model, this region has a radius of ~ 6 pc, while for 
the STARS model the radius is only ~ 2 pc, the 
shell is much narrower and the densities are higher. 
This is because the RSG phase is much shorter in 
the STARS models, even though roughly the same 
amount of stellar mass is lost as in the MM2003 
models. The ram pressure of the RSG wind is less 
than the thermal pressure of the hot, shocked bub- 
ble produced during the main-sequence stage. As 
a result, there is some back filling as the hot dif- 
fuse gas pushes back towards the star. The lack of 
ionizing photons at this stage (see Fig. [3]) , and en- 
hanced opacity due to the RSG wind and neutral 
shell, leads to recombinations in the outer H II re- 
gion at radius R ~ 27 pc. This is most pronounced 
for the MM2003 model, since the inner absorbing 
RSG shell contains more neutral atoms than the 
shell formed in the STARS case. When thermal 
conduction is included in the simulations, evap- 
oration of cold material from both the RSG shell 
and the outer H II region enhances the density and 
thus the cooling rate in the hot shocked wind bub- 
ble. In the MM2003 case, this results in a much 
smaller hot bubble (~ 19 pc as opposed to ~ 27 pc 
in the simulation without conduction). The cir- 
cumstellar medium in the STARS simulations is 
less affected by thermal conduction, possibly due 
to the shorter timescales in the RSG phase. 

In Figure [7] we show the corresponding circum- 
stellar structures for the 60 M© models at the end 
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Table 1 

Radii and masses of the outer H I shell at the end of MS phase 



Mass 


Model 




M ism 


M 




pc 


M Q 


40 


MM2003 N 


24.6 


2.0 x 10 5 


40 


MM2003 R 


27.5 


2.8 x 10 5 


40 


STARS 


29.0 


3.3 x 10 5 


60 


MM2003 N 


25.6 


2.2 x 10 5 


60 


MM2003 R 


29.4 


3.4 x 10 5 


60 


STARS 


29.4 


3.4 x 10 5 




Radius [pc] 



5 10 
Radius [pc] 



10 

12 2 

14 -5 



Fig. 5. — Density and temperature at 10 6 yrs of evolution for the 4OM model from STARS. Left-hand 
panels: model without thermal conduction, right-hand panels: model with thermal conduction. The density 
plots show total density (solid line) and ionized density (dashed line). The temperature plots also show the 
internal energy e = p/[y — 1] (dashed line). 
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Fig. 6. — Density and temperature at the end of the RSG phase for the 40 M & models. Left-hand panels: 
MM2003 models, right-hand panels: STARS models. Top row: models without thermal conduction, bottom 
row: models with thermal conduction. The density plots show total density (solid line) and ionized density 
(dashed line). The temperature plots also show the internal energy, e = p/[y — 1] (dashed line). 
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Fig. 7. — Same as Fig. [6] but for the 60 M & models. 
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of the RSG phase. The MM2003 and STARS 
results look quite different. This is because the 
MM2003 model experiences several episodes of en- 
hanced mass loss, whereas the STARS model un- 
dergoes a single intense outburst. The ram pres- 
sure of the slow wind in the MM2003 case fluctu- 
ates and this leads to acoustic waves propagating 
from the inner edge of the hot bubble outwards. 
In Figure [7] the waves have yet to reach the outer 
edge of the bubble, but once they do, they will 
reflect from the dense neutral shell and pass back 
through the bubble and reflect backwards and for- 
wards through the hot shocked gas, forming lo- 
calized heated and compressed regions where op- 
positely directed waves shock against each other. 
There is no clearly defined shell of RSG material 
in the 60 M Q models, because the duration of the 
most intense period of mass loss is very short and 
just before the onset of the WR phase (see Fig.Q]), 
hence the dense material has not had time to move 
far away from the star before the fast wind starts. 

We are particularly interested in the radii of 
the dense, neutral RSG shells formed in these 
models. These shells will be impacted and bro- 
ken up by the fast WR winds and so have im- 
portant consequences for the observables (X-ray 
and optical emission) of the WR wind-blown bub- 
bles. The shells are most evident in the 40 M 
models, reaching radii of ~ 6 pc and ~ 2 pc for 
the MM2003 and STARS models, respectively, by 
the end of the RSG phase. These distances re- 
flect the duration of the RSG stage in each case: 
5.6 x 10 5 yrs for the MM2003 model as opposed 
to 3.1 x 10 5 yrs for the STARS model. Also, the 
wind velocities in the STARS model are allowed 
to fall below 30 km s _1 and in general are lower 
that those of the MM2003 models (see Fig.©. We 
have not presented results for the MM2003 mod- 
els with stellar rotation because they are qualita- 
tively similar to the results for the other models 
and result in a dense neutral shell being formed at 
~ 2 pc, since the RSG phase in this model only 
lasts 1.0 x 10 5 yrs. We continue the simulations in 
ID until the end of the WR phase but we do not 
show these models because we are interested in the 
formation of instabilities during the fast wind-slow 
wind interaction, which do not occur in ID. 



4.2. 2D Results 

At the end of the RSG phase, we remap our 
ID, spherically symmetric distributions of den- 
sity, momentum and energy onto a 2D axisymmet- 
ric grid using a volume-weighted averaging pro- 
cedure to ensure conservation of these quantities. 
The 2D grid is, of computational necessity, coarser 
than the ID grid, and we also zoom in on the 
inner ~ 10 pc of the simulation, since this is 
where the hydrodynamic interaction between the 
fast WR wind and the slow RSG wind will take 
place. Also, the observed WR bubbles all have 
radii 10 pc or less. We continue to follow the 
radiation-hydrodynamic evolution of the circum- 
stellar medium using the stellar wind parameters 
and the tabulated ionizing photon rate as a func- 
tion of time shown in Figures. [TJ [2j and[H For all 
2D simulations, we present only one quadrant of 
the calculation, since the other is the same due to 
symmetry. 

4.2.1. 40 M models 

In Figure [5] we present grayscale images for 
the wind-wind interaction for the 40 M STARS 
model without thermal conduction at four differ- 
ent times corresponding to when the interaction 
region has reached distances 2, 4, 6, and 8 pc 
from the central star. These times are, respec- 
tively, 39800,65900,86000, and 104000 yrs after 
the onset of the WR wind. In this figure, we see 
the disruption of the RSG shell by the fast wind 
and the formation of clumps. The RSG material 
consisted initially of a thin, dense shell of partially 
ionized material located ~ 2 pc from the star. The 
fast WR wind accelerates down the density gradi- 
ent and slams into the dense shell. The WR wind 
is shocked and heated by this interaction while the 
dense shell is shocked and compressed. Vishniac 
instabilities of the linear thin-shell variety break 
the shell into de nse clumps, whi c h then travel ra- 
dially outwards ( Vishniac 19831 : Vishniac fc Rvu 
19891 ) . The clumps have a cometary form, with the 



dense, cold (10 4 K photoionized) head pointing to- 
wards the central star and tails extending several 
parsecs radially outwards, showing that the heads 
are traveling more slowly than the diffuse material 
that surrounds them. Interactions between clump 
material and the shocked fast wind create a pat- 
tern of multiple interacting shocks, which heat the 
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Fig. 8. — Wind-wind interaction for the 40 Mq STARS model without thermal conduction. From left to 
right the evolution times are 39800, 65900, 86000, and 104000 yrs after the onset of the WR wind. Top panels 
show ionized number density and bottom panels show temperature, both on logarithmic scale. 
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Fig. 9. — Same as Figure [8] but for the 40 M & MM2003 model without stellar rotation and without thermal 
conduction. From left to right the evolution times are 10600, 22600, 32650, and 36600 yrs after the onset of 
the WR wind. 
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Fig. 10. — Same as Figure M but for the 40 M© MM2003 model with stellar rotation without thermal 
conduction. From left to right the evolution times are 6800, 8800, 12800, and 16900 yrs after the onset of the 
WR wind. 
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ablated clump gas to X-ray temperatures. For this 
particular case the RSG shell breaks up at ~ 5 pc 
(see Fig.!]). 

Figure [5] shows the results obtained when the 
stellar wind parameters come from the MM2003 
40 Mq model without rotation and without con- 
duction. The wind-wind interaction region is 
again shown at 2, 4, 6 and 8 pc from the central 
star, corresponding to 10600, 22600, 32650, and 
36600 yrs after the onset of the WR wind. For 
this model, the RSG dense wind material, trav- 
eling at ~ 30 km s _1 , is in the form of a broad 
shell about 6 pc from the star at the start of the 
WR stage. The interaction of the WR wind with 
this shell occurs only a few times 10 3 yrs after the 
onset of the fast wind. Compared to the STARS 
model, the MM2003 models experience Rayleigh- 
Taylor instabilities, which disrupt the accelerating 
WR shell before it interacts with the main shell of 
RSG material, forming small clumps at < 2 pc. 
The main interaction of the WR shell with the 
RSG shell occurs at 6 pc and compresses the RSG 
material into a very thin, dense shell, which then 
corrugates and starts to break up at around 8 pc 
due to the linear thin-shell instability. 

Figure [10] shows the results from MM2003 
model with stellar rotation but without conduc- 
tion at 6800, 8800, 12800, and 16900 yrs after the 
onset of the WR wind. In this model the interac- 
tion takes place much further from the star and 
the shell of RSG material is not so dense, so the 
interaction does not produce such violent insta- 
bilities and the shell remains essentially intact, 
being swept up and compressed by the fast wind 
but not suffering the complete disruption seen in 
the STARS case discussed above. For this model, 
the intense mass-loss phase has a relatively short 
duration, ~ 10 5 yrs, and there are several mass- 
loss episodes before and after the main ejection 
event. Rather than a RSG stage, the extremely 
strong mass loss (M w ~ 10 _36 Mq yr _1 ), together 
with the fact that the stellar effective temperature 
is relatively high (T e g ~ 10 4 K) for part of this 
period, indicate a Luminous Blue Variable stage, 
with stellar wind velocities V w > 200 km s _1 . The 
densest wind material travels quite a large dis- 
tance from the star before the onset of the WR 
stage and so the interaction with the fast wind, 
when it occurs, does not lead to complete breakup 
of the shell. 



We find that thermal conduction has little dis- 
cernible effect either on the density or the tem- 
perature structure of the 2D simulations and the 
results are essentially the same as those shown in 
Figures [5] M and [TDJ In Figure Qj] we present one 
set of results including thermal conduction to il- 
lustrate this. In this figure we show the results 
of the 40 M Q MM2003 simulations (without rota- 
tion) at 36600 yrs after the onset of the WR wind. 
There are slight differences between the two sets of 
results, but the short timescales of the wind- wind 
interaction stage compared to cooling timescales 
in the 10 6 K gas mean that there is no major effect 
on the hydrodynamics due to thermal conduction 
during this stage. 

We have run further set of STARS simula- 
tions to investigate the effects of the radius of the 
wind injection zone and the grid resolution on the 
number and distribution of clumps formed during 
shell disruption. Increasing the grid resolution to 
800x1600 cells while maintaining the same num- 
ber of cells in the wind injection region (radius 30 
cells) produced the same results as those presented 
here. Decreasing the size of the wind injection re- 
gion (radius 10 cells) did produce fewer clumps but 
increasing the radius of the wind injection region 
to 50 cells gave the same results as in Figure [8] 
The clumps formed in our simulations have diam- 
eters between 15 and 50 cells and are photoionized, 
hence are fully resolved. 

4.2.2. 6OM models 

In Figure [12] we show the corresponding results 
for the 60 Mq STARS model without conduction. 
Again, we show the wind-wind interaction at 2, 
4, 6, and 8 pc, corresponding to evolution times 
of 27150, 43200, 57200, and 71250 yrs after the on- 
set of the WR wind. The dense RSG (or LBV) 
material is swept up into a thin shell by the fast 
WR wind. Linear thin-shell instabilities are in- 
duced in the swept-up material and break the shell 
into many cold (< 10 4 K), dense (> 1 cm~ 3 ) 
clumps which drift outwards, embedded in the 
hot, shocked WR wind. The morphologies of the 
40 Mq and 60 Mq STARS models are very simil- 
iar. 

For the MM2003 60 M model without rota- 
tion, there are several episodes of enhanced mass 
loss during the RSG/LBV stage and even in the 
WR stage there are short-lived outbursts (see 
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Fig. 11. — Ionized density (Top panels) and temperature (bottom panels) during the interaction of the fast 
WR wind with the slow, dense RSG wind for the 40 Mq MM2003 model without stellar rotation. Left panel: 
model without thermal conduction Right panel: with thermal conduction. Both numerical results are at 
36600 yrs after the the onset of the WR stage. 
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Fig. 12. — Wind-wind interaction for the 60 M STARS model without thermal conduction. From left to 
right the evolution times are 27150, 43200, 57200, and 71250 yrs after the onset of the WR wind. 
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Fig. 13. — Same as Figure [T2l but for the 60 M & MM2003 model without stellar rotation and without thermal 
conduction. From left to right the evolution times are 6800, 14800, 24800, and 30900 yrs after the onset of 
the WR wind. 
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Fig. 14. — Same as Figure fT2l but for the 60 M Q MM2003 model with stellar rotation and without thermal 
conduction. From left to right the evolution times are 19450, 21450, 23450, and 25500 yrs after the onset of 
the WR wind. 



22 



Fig. [T]). This leads to the formation of succes- 
sive shells of material traveling with different ve- 
locities away from the star shown in Figure 1131 
The most intense mass-loss episode shown in Fig- 
ure Q] is followed by a reduction in mass-loss rate 
and increase in wind velocity (see Fig. [2]), which 
sweeps up, compresses and accelerates the ejected 
material. This leads to the formation of Rayleigh- 
Taylor instabilities, which disrupt the shell once it 
has reached ~ 10 pc from the star. The clumpy 
shell is traveling quite quickly and soon leaves 
the region of interest (i? < 10 pc). Other, less 
dense, shells follow, which suffer instabilities but 
do not contain enough material to form dense, cold 
clumps. In Figure [T3] we show the wind- wind in- 
teraction and the formation of different shells and 
shock patterns as a result of the variability of the 
mass- loss and wind velocity up to 30,900 yrs after 
the onset of the WR wind. 

For the MM2003 60 M Q model with rotation, 
from Figure [1] we note that there is no period 
of particularly enhanced mass loss. In this case, 
therefore, although the RSG is swept up by the 
faster WR wind, the swept-up material is not par- 
ticularly dense and does not form a thin, dense 
shell nor any instabilities or clumps. In Figure [T4l 
we show the density and temperature in the evolv- 
ing CSM up to the point when the swept-up ma- 
terial has reached a radius of 8 pc from the star, 
which occurs 25,500 yrs after the onset of the WR 
stage. 

4.3. Comparisons with previous works 

Differences between our results and the work 
of previous authors can be attributed to several 
factors: the stellar evolution models, the numeri- 
cal method, and the initial conditions and physics 
that is included in the numerical simulations. In 
this section we discuss each of these aspects in 
turn. 

Much of the previous numerical work on the 
evolving circumstellar media around massive stars 

s tellar e vo- 
Garcia-Segura et al.l (Il996allbh . 



has used the sam e 35 M w and 60 M P) 
lution models as 



which are those developed bv lLanger et al.l ([1994D . 
These models were used in their full t im e-dependent 



form b y Garcia-Segura et all dl996alrbT) 
(|2003l ' l2006l ). and iDwarkadasI (|2007). 



Frever et al 



For both 

stellar masses considered, the stars evolve to be 
red supergiants (T c g < 4800 K). One difference 



between these models and the ones we use in 
this paper is the stellar w ind velocity. In the 



Garcia- Segura et al. (|l996al lbh models, the stellar 



wind velocity remains high, above 2500 km s , 
throughout the main-sequence stage, whereas in 
all of our models we see a significant decrease 
from an initial ~ 3000 km s _1 to ~ 1000 km s _1 . 
This affects both the pressure and size of the 
hot, shocked wind bubble at the end of the main- 
sequence stage. In the case of their 60 M Q model, 
there is a "H-rich WN" stage immediately af- 
ter the main-sequence stage, in which more than 
10 M© of material is lost by the star at high ve- 
locity over the space of a million years. None of 
the models we have studied have a similar evo- 
lutionary stage. This mass is moving too fast to 
remain close to the star during the subsequent evo- 
lutionary stages, and thus does not contribute to 
structure formation in the circumstellar medium. 
Other authors have used the STARS stellar evo- 
lution models ( Eldridge et al. 20061) . where the 



only differences between that work and ours is 
that we have corrected the stellar wind veloci- 
ties (upw ards) by missing num erical factors, and 



those of Schallc r et al 
60 M Q (|van Marie et al 



19921) for 40 M Q and 
20051 120071) . which arc 



very similar to the MM2003 model s with o ut ro - 
tation. However, Ivan Marie et al.l (120051 . 120071 ) 



approximates the evolution with a three-constant- 
winds scenario, and so does not follow the episodic 
ejections of material which are a main character- 
istic of the post-main-sequenc e evolution of star s 
with initial masses > 60 M© (|Humphrevsll2010h . 
By following the full time variation of the stellar 
wind parameters, we are able to study the for- 
mation and interaction of multiple shells in the 
circumstellar medium in the LBV and WR stages. 

All of the numerical simulations by previous 
authors use well-tested hydrodynamical schemes 
and so we do not expect major differences to 
arise from differenc es between codes. The cod es 
used include ZEUS ( Garcia-Segura et al. 1996alfb ; 



Ramirez-Ruiz et al.l 2005; van Marie et al. [2005 



Eldr idge et al.ll2006 : van Marie et al. l2007tlChita et al.l 
l2008tlPerez-Rend6n et al.ll2009h. VH1 (IDwarkadasI 



l2007h and a code due to 



and lYorke fc Weld (|l996l ) (|Frever et al.l |200 



Yorke fc Kaisigj (1199 



2006). Some authors choose spherical coor- 
dinates for their two-dimensional models (e.g. 



Garcia-Segura et al.lll996allb1 ; Ivan Marie et al.l2005l 
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Dwarkadad l2007t Ivan Marie et al.l 120071 ). while 
others us e cylindrical r — z coordinates, 



as we 



do (e.g., IFrever et all 120031 . l200d ). Each coordi 



nate system has its advantages and disadvantages. 
In 2D spherical coordinates, the cell sizes increase 
with increasing radius, while in r — z coordinates 
it is more difficult to represent a spherical wind- 
injection zone. 

One key ingredient, which can strongly affect 
the results of the numerical simulations, is the 
inclusion (or not) of radiative transfer and the 
modeling of the H II region that will develop 
around the star during the main-sequence and 
WR stages. Of the p revious nume r ical work w e 
have mentioned , only I Frever et al" ( 2003, [2 006), 



van 



Marie et al. (|2005l . 120071 ) and IChita et al 



20081) (for &Y2 M e star) include the effect of ion- 
izing photons o n the evolution of the ci r cums tellar 
medium. Both van Marie et al.l (|2005l 120071 ) and 
Chita et al. ( 20081) use a simple Stromgren radius 
approach to the radiative transfer, whereby gas 
is either ionized or neutral depending on whether 
the ionizing photons reach it or not (it is not clear 
whether these authors take c ollisional io nization 
into account). IFrever et al. (l2003l 120061) uses a 
ray-tracing method and includes collisional ioniza- 
tion. Our radiative transfer method (see § 12 - 1[) al- 
lows us to follow gas at intermediate stages of ion- 
ization. For example, once the RSG stage begins, 
the ionizing photons do not reach the outer radii 
of the nebula, but the recombination timescale 
is of order 10 5 /n yrs and thus the gas does not 
recombine immediately. This has important con- 
sequences for the pressure across the bubble. In 
our results, even once the ionizing photons have 
switched off, we find roughly constant pressure 
right across the bubble, from the inner edge of the 
hot, shocked wind through the recombining H II 
region and the outer neutral shell, since acous- 
tic waves smooth out the pressure distribution as 
the H II region slowl y recombines. Howeve r, the 
results pre sented by | van Marie et al. (I2004L Fig- 
ure 4) and van Marie et al.l ( 2005 . Figure 1) show 
a drop in pressure in the ex-H II region where the 
number of particles has suddenly been reduced by 
half due to instantaneous recombination. This will 
lead to the formation of spurious shells of material 
as shock waves propagate into the depressurized 
region from both sides. 

The expansion of a stellar wind bubble into an 



H II region cannot simply be modeled by a bub 
ble expanding into a unifo r m medium at 10 4 K 
(IGarcfa-Segura et al.lll996aT lb1: lRamirez-Ruiz et al 



2005), since the density of the photoionized gas in 



the H II region becomes lower with time as the 
expansion proceeds. That is, not only is the stel- 
lar wind bubble expanding, but the H II region is 
expanding as well. 

In our 2D models we see the Rayleigh- Taylor 
and linea r thin-shell instabilities rep or ted by other 
authors (IGarcfa-Segura et al.ll996al lbl; IFrever et al. 
2003Ll2006UDwarkadasl2007t ) for the wind- wind in- 
teraction stage. The details very much depend on 
the stellar evolution model being used. Because 
we do not model the main-sequence stage in 2D we 
do not see the H II shadowin g ins tability reported 
bv IFrever et all (|2003l liiooi ) and lArthur fc Hoard 
(|2006l ). 

We remark that the initial uniform ambient 
medium density of our models is no = 100 cm" 3 , 
which is higher than that used by all the previous 
authors we have mentioned. The justification for 
our value is that a massive star will live a large 
fraction of its lifetime in the vicinity of the molec- 
ular cloud where it was born. As the H II region 
around the star evolves, the density of the pho- 
toionized gas falls with time. For example, elec- 
tron densities in the Orion nebula (age about one 
millio n years) range from 100 cm" 3 dFelli et al. 
19931) to > 1000 cm" 3 (jOsterbrock fc Flather 



1959), while those in the Carina nebula (which 



harbors several high-mass main-sequence stars 
and post-main-sequence st ars) vary between ~ 
30 cur 3 and ~ 350 cm" 3 (|Mizutani et al.ll2002t ). 
By the end of the main sequence, the density in 
our photoionized gas has fallen to < 10 cm" 3 . The 
stellar wind bubble evolves within this changing 
environment, with the position of the inner wind 
shock regulated by the pressure in the H II re- 
gion. Models in which the external med i um is 
initially less dense (e.g., Frever et al. 2003, 2006; 



van Marie et al.ll2005l 120071 ) produce larger, more 



diffuse H II regions and consequently larger stel- 
lar wind bubbles than in our models. Simula- 
tions with a lower density medium and without 
radiative transfer produce large stellar wind bub- 
bles with thin ou ter shells of swept-up material 
(|Dwarkadasll2007l ) instead of broad shells of pho- 
toionized gas. 
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4.4. Energy evolution 

In Figure [15] we show the evolution of neutral 
and ionized gas kinetic and thermal energy frac- 
tions with respect to the total integrated stellar 
wind mechanical energy, Ey,(t) = J M w V%/2dt, 
resulting from the numerical simulations from the 
40 M Q STARS models with and without thermal 
conduction. We do not include the energy of 
the stellar ionizing photons in this calculation be- 
cause nearly all of the interstellar UV energy will 
be radiated away in forb idden-line processes (see 
Dvson fc Williams! 1199/j Chapter 7). At the be- 
ginning of the main-sequence stage, the energy 
fraction profiles show some variations due to the 
formation and e xpansion of the in i tial H II regi on 
around the star (jFrever et al.ll2003l:lArthuH[20 07'FI. 
For most of the main-sequence stage, the energy 
fractions of both kinetic and thermal ener gy are 
roughly constant (jFrever et al.ll2003l [20061) . with 
the thermal energy being lower in the model with 
conduction due to enhanced cooling in the hot 
bubble because the material evaporated from the 
dense outer shell increases the density, and hence 
the cooling rate, in the outer part of the bubble. 
It is evident that the outer neutral shell dominates 
the kinetic energy fraction of these models, while 
the hot, shocked bubble dominates the thermal 
energy. 

At the end of the main-sequence stage, the neu- 
tral shell remains coasting along but the hot bub- 
ble suffers a drastic drop in pressure when the 
(subsonic with respect to the hot bubble) RSG 
wind starts. This can clearly be seen in Figure [T5l 
where, starting at about 4 x 10 6 yrs, the thermal 
energy fraction dips while the kinetic energy frac- 
tion remains roughly constant. With the onset of 
the fast WR wind, the thermal energy increases 
enormously once the high velocity gas shocks and 
thermalizes into a new pressure-driven hot bubble. 
In contrast, the kinetic energy becomes shared be- 
tween the ionized and neutral components. This is 
due to photoionization of part of the outer neutral 
shell and the repressurized hot bubble rejuvenat- 
ing the expansion of the outer swept-up material. 

Between them, the thermal and kinetic energy 



fractions add up to substantially less than 1, par- 
ticularly during the main-sequence stage. The re- 
maining energy is lost as radiation both from the 
hot bubble and from the H II region and neutral 
shell which surround the wind bubble. 

5. Comparision with observations 
5.1. Optical observations 

Our simulated WR wind bubbles bear morpho- 
logical similarity to optical images o f observed WR 



nebulae such as RCW 58 and S 308 (Grucndl et al 



2000). Instabilities in the swept-up RSG material 
shell lead to the formation of dense, cold clumps, 
which can be compared to t he Ha knots, for exam - 



Gruendl et aD (|2000h 



3 Note that iFrever et all d2003l. 1200611 do include the stellar 
UV energy in their calculations, which leads to lower gas 
thermal and kinetic energy fractions by more than an order 
of magnitude. 



pie, in images of RCW 58. 
divide a sample of WR bubbles into morphological 
types depending on the offset of the Ha and [O III] 
emission. Type I has no measurable offset between 
the Ha and [O III] fronts, Type II morphology has 
an Ha front trailing close behind an [O III] front of 
similar shape, while for Type III the Ha trails far 
behind a faint [O III] front and can differ appre- 
ciably in shape. Type IV morphology has a faint 
[O III] front with no Ha counterpart. Under this 
classification system, RCW 58 is Type III, while 
S 308 is Type II. While the Ha is a recombination 
line produced in gas with temperatures ~ 10 K, 
[O III] comes from hotter gas > 10 4 K due to 
either photoionization by harder UV photons or 
collisional excitation behind shock waves. In our 
simulations, the dense clumps that result from the 
instabilities which break up the swept-up shell are 
cold (T < 10 4 K) and are photoionized on the side 
closest to the star. These clumps will therefore 
be sources of Ha emission and, moreover, lag be- 
hind the remainder of the shocked shell, which has 
temperatures T > 10 4 K and could be a source of 
[O III] emission. The images presented in section 
4.21 were chosen to facilitate comparison with the 
ring nebula S308, which has an optical radius of 
- 8 pc. 

The initial abundances in the stellar evolu- 
tion models from MM200 3 and STARS are Solar 
(|Anders fc Grevessd[l989T ) 3ut towards the end of 
the RSG phase the stellar surface abundances be- 
gin to change, becoming enriched with nitrogen in 
particular. This surface material is expelled from 
the star in the form of stellar winds during the 
period of enhanced mass loss and can be expected 
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Time [Myr] Time [Myr] 

Fig. 15. — Kinetic and thermal energy as fractions of the time-integrated wind mechanical energy as a 
function of time for the 40 M Q STARS models. Left panels: Gas kinetic energy as a fraction of total wind 
mechanical energy. Right panels: Gas thermal energy. Solid thick lines: Total gas kinetic (thermal) energy. 
Solid thin lines: Neutral gas kinetic (thermal) energy. Dotted lines: Ionized gas kinetic (thermal) energy. 
Top panels: simulations without thermal conduction, bottom panels: simulations with thermal conduction. 
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Table 2 

Stellar and nebular parameters from selected WR wind-blown bubbles 



WR bubble 


Distance a 


Radius' 3 


logM c 


V c 


log N/O d 




ye 






kpc 


arcmin 


M yr- 1 


km s -1 




10 21 cm" 2 


K 


10 34 erg s _1 


S308 


1.5±0.3 


19.4 


-4.12 


1720 


0.22 


1.1 


1.1 x 10 6 


1.2 


NGC 6888 


1.8±0.5 


8.6 


-4.02 


1605 


0.27 


3.1 


1.3 x 10 6 


3 


RCW58 f 


3.0 


4.9 




975 


-0.30 









a Distances were obt ained from : Chu et all (J2003) 
WR136 (NGC 6888), andfch] (Il982h - WR40 (RCW58). 



WR6 (S308), lAbbott et all (119861) 



1 Gruendl et al. ( 


2000) 


'Prinia et al. ( 


1990|). 


c Esteban et al. ( 


1992) 



e Obtained from soft X-ray observations ( Chu et al. 2003 ; Wrigge et al. 20051) . 
f RCW58 has not been detected in X-ray observations ( Gosset et al.ll2005j ) . 



to enrich the circumstellar nebula. Indeed, spec- 
troscopic abu ndance determination s of Wolf-Rayet 
ring nebulae (jEsteban et al.lll992n show that the 
nitrogen-to-oxygen ratio, log(N/0), can even be- 
come positive in this circumstellar medium, while 
the Solar value is log(N/0) = -0.9 (see Table©. 

In Figure [16] we plot the nitrogen-to-oxygen 
abundance ratio for stellar material expelled since 
the onset of enhanced mass loss for the different 
stellar evolution models used in this paper, us- 
ing the model stellar surface abundances as being 
representative of the stellar wind chemical com- 
position. These are mass-weighted average val- 
ues over this time period and do not include ma- 
terial expelled during the main-sequence stage. 
The main-sequence stellar wind material is spread 
throughout an extended, diffuse bubble of radius 
up to some tens of parsecs, while the RSG and 
Wolf-Rayet wind material will form a circumstel- 
lar nebula around the Wolf-Rayet star. We com- 
pare our calculated abundance ratios to those de- 
termined for the three Wolf-Raye t ring nebulae 
S308, NGC 6888 and RCW58 bv lEsteban elaT 



(|l992j ) and listed in Table [2] It is clear from these 
figures that rotation significantly affects the abun- 
dances at the stellar surface, and hence in the stel- 
lar wind, throughout the main-sequence stage so 



that by the start of the enhanced mass-loss pe- 
riod the nitrogen-to-oxygen ratio is already posi- 
tive. For the evolution models without rotation, 
the abundances at the start of the RSG stage are 
still Solar, since there has been no mixing with 
enriched material from the stellar core up to this 
point. 

From Figure [16] we see that the abundance ratio 
for RCW 58 can only come from material expelled 
in the early stages of enhanced mass loss. The 
S 308 and NGC 6888 values could be reproduced 
by various models, except that of the MM2003 
40 M Q without rotation. For the STARS models, 
however, these abundances are achieved only once 
the Wolf-Rayet stage has been entered and WR 
wind material has been added to the circumstellar 
nebula. For the MM2003 models with rotation, 
the nitrogen-to-oxygen abundance ratio in the cir- 
cumstellar medium is predicted to be far higher 
than the observed values throughout the RSG and 
WR stages. These findings are broadly consis- 
tent with the classification of the central star of 
RCW 58 as WNL, while the central stars of S 308 
and NGC 6888 are WNE. 
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Fig. 16. — Nitrogen-to-oxygen abundance ratio 
in the circumstellar nebula as a function of time 
since the end of the main-sequence stage for the 
stellar evolutionary models used in this work. 
Top panel: 40 Mq models, bottom panel: 60 Mq 
models. In each panel the solid line is for the 
STARS stellar evolution models, the dashed line 
is for the MM2003 models without rotation, while 
the dotted line is for the MM2003 models with 
stellar rotation. Also shown are thin horizontal 
lines corresponding to the spectroscopically deter- 
mined values for the three Wolf-Ray et ring nebulae 
S308, NGC 6888 and RCW58 from|Esteban etal 



(|1992I ). The triangles on each line mark the point 
where each stellar model enters the Wolf-Rayet 
phase. 
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Fig. 17. — Differential emission measure for the 
models presented in Figs. [51 |H and [10] when the 
wind-wind interaction is at 8 pc. Top panel: 
models without conduction, bottom panel: models 
with conduction at similar points in their evolu- 
tion. In each panel, dashed line — STARS model, 
solid line: MM2003 without stellar rotation, dot- 
ted line: MM2003 with rotation. The temperature 
bin width is 0.04 dex in each case. 
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5.2. X-ray observations 

There are around 200 WR stars in our Galaxy, 
and only two of them have nebulae detected in 
diffuse X-rays: S308 and NGC6888 (iBochkarev 



1988 [ IWrigge et al.lll994t IWriggel ll999t IChu et al 
2003'; IZhekov fc Parkl l201ll) . The X-ray emis- 



sion is very soft and the diffuse emission from 
S 308 can be fit with a one temperature compo- 
nent model with (T = 1.1 x 10 6 K), enhanced 
nitrogen abundance and X-ray luminosity Lx — 
1.2 x 10 34 erg s _1 in the 0.25-1.5 k eV energy 



band , assuming a distance of 1.8 kpc (jChu et al 



20031 ). NGC6888 ASCA and ROSAT X-ray ob- 
servations are best fit by a two-component model 
with the dominant component at 1.3 x 10 6 K and 
the weaker component at 8 x 10 6 K and a to- 
tal X-ray luminosity Lx — 3 x 10 34 erg s _1 in 
the energy rang e 0.4-2.4 keV, assuming adis- 
tance of 1.8 kpc (jWrigge et al.lll994 Il998i 120051) . 



The recent SUZAKU observations of NGC 6888 
(|Zhekov fc Parkl201ll) confirm that ninety percent 
of the observed X-ray emission is in the soft 0.3- 
1.5 keV energy range. In both bubbles, the X-ray 
emission is internal to the optical [O III] emission. 
Such soft X-ray emission indicates that the obser- 
vations will be strongly affected by the absorbing 
column density, ATh, along the line of sight. This 
could be a reason for the non-detection of other 
wind-blown WR bubbles such as RCW58, since 
these objects are found in the plane of the Galaxy 
where Nh will be greatest. The two bubbles with 
detected X-ray emission are relatively close by 
whereas WR40, the central star of RCW58, is 
more than twice as distant (see Table [2]). 

5.2.1. X-ray simulations 

From our numerical simulations, both in one 
and two dimensions, we can calculate the expected 
X-ray luminosity for any energy range at any time 
in the simulation. To do this efficiently, we use the 
temperature and density in each cell of the hydro- 
dynamic simulation to construct the array of Dif- 
ferential Emission Measures (DEM), i.e. the emis- 
sion measure (EM. = J n e riidV) that corresponds 
to each bin of a preassigned temperature array for 
a given numerical simulation. Note that we mask 
out the unshocked thermal wind in this procedure. 
The total X-ray spectrum is then found by sum- 
ming the spectra of each temperature bin and the 



X-ray luminosity is then found by summing over 
the required energy range. This method is much 
more efficient than summing the spectrum (or lu- 
minosity) of each individual hydrodynamic cell, 
since we find that 100 temperature bins are per- 
fectly adequate, while the numerical simulations 
consist of ~ 6000 cells for the one-dimensional 
models and > 10 5 cells for the two-dimensional 
models. In Figure [17] we show the corresponding 
DEM histograms for the 40 Mq results presented 
in Figures. [Hini and 1 101 when the wind- wind inter- 
action is at 8 pc, together with their counterparts 
that include thermal conduction. 

The first thing to notice about Figure [17] is that 
most of the mass has temperatures T < 10 6 K. 
This means that the spectra will be very soft and 
the absorbing column of neutral hydrogen towards 
such objects will play a very important role in 
what is observed. We see no systematic difference 
between models without thermal conduction and 
their equivalent models with conduction beyond 
the differences that can be attributed to slightly 
different volumes of emitting gas. The main dif- 
ference we see is between models with complete, 
if distorted, shells (those of MM2003) and those 
where the shell has already broken up into clumps 
(those of STARS). The models with clumps have 
a much smaller differential emission measure for 
the hot gas (T > 10 5 5 K) than those where the 
whole shell is being shocked, because the hot gas 
in these cases is more diffuse and can flow around 
the clumps and out of the computational domain. 

Once we have the DEM arrays, the X-ray 
spectra are calculated assuming collisional ion- 
iza tion equilibrium (CIE) and the ion fractions 



of iMazzotta et al.l (]1998l) . The line spectrum 



uses the APED database of CIE line intensities 



([Smith et al.ll200ll ). while the continuum is calcu- 
lated with our own code, taking into account free- 
free, free-bound and two-photon emission. The 
chemical abundances can be varied according to 
the requirements of the models. 

Typical electron densities in the X-ray emitting 
gas for our simulations are n e ~ 0.1 cm~ 3 and 
temperatures are in the range 1-4 x 10 6 K (see 
Fig. [5] through Q3] and [T7]) . For these parame- 
ters, we can assess the closeness of the different 
elements (C, N, O, N e, etc.) to ionization equ ilib- 
rium using Figure 1 of lSmith fc Hugh es (2010). At 
10 6 K, for nitrogen (the main X-ray spectral char- 
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acteristic of S 308, Chu et al . 2003) to be in ioniza- 



tion equilibrium requires n e t > 8 x 10 cm s , 
whereas at 4 x 10 6 K, ionization equilibrium for 
nitrogen requires n e t > 1 x 10 11 cm -3 s . For 
our typical electron density, these correspond to 
timescales of 253,000 yrs and 32,000 yrs respec- 
tively. Thus, it could be argued that a non- 
equilibrium ionization treatment is required for 
the X-rays in the shocked stellar wind interaction 
region, but this is beyond the scope of this paper. 

In Figures [18] and [19] we show the X-ray lumi- 
nosity as a function of time in the WR phase for 
the 40M(D models, assum ing Solar abundances 
(|Anders fc Grevessd Il989h . Figure HH plots the 
soft-band (0.25-1.5 keV) and hard-band (1.5- 
7 keV) X-ray emission for the ID STARS models 
with and without thermal conduction, together 
with the STARS 2D model emission with thermal 
conduction. From this figure we see that the ID 
models produce more intense X-ray emission, by 
an order of magnitude, than the 2D model. This 
can be understood by referring to Figure |8] where 
we see that the RSG shell breaks up into clumps 
and the low-density, hot, shocked WR wind flows 
around them. The X-ray emission starts to dimin- 
ish after ~ 70,000 yrs because the hot gas flows 
off the 10 pc 2D grid. 

In the ID models, the WR wind is confined by 
the shell of RSG material, which cannot break up 
because of the geometrical restriction. The soft 
X-ray emission of the two ID models is very sim- 
ilar until 1.5 x 10 5 yrs, when the thermal conduc- 
tion at the edge of the hot bubble starts to play 
a dominant role, raising the density and lowering 
the temperature of the hot gas here and thereby 
increasing the contribution of the soft X-rays. In 
the ID model without conduction, the gradual in- 
crease in the stellar wind velocity raises the tem- 
perature of the hot, shocked gas and finally re- 
sults in the hard X-rays being the dominant com- 
ponent. Note that, for the ID models, we are sum- 
ming the X-ray contribution from the whole com- 
putational grid (except the thermal wind injection 
region): the contributions from the diffuse, hot, 
shocked main-sequence wind to the soft and hard 
X-rays are, respectively, 2.1 x 10 31 erg s -1 and 
2.9 x 10 31 erg s -1 for the model without conduc- 
tion and 1.6 xlO 31 ergs -1 and 2.1 xlO 31 ergs -1 for 
the model with thermal conduction (see Fig. [6j. 

Figure [19] shows the complete time evolution 



E 
s 




Time [10 = yr] 

Fig. 18. — Variation of X-ray luminosity with time 
in the WR phase for the STARS 40M Q model 
for ID and 2D simulations, assuming Solar abun- 
dances. Solid line: ID simulation with thermal 
conduction, dashed line: ID simulation without 
conduction, dotted line: 2D simulation without 
conduction. In each case the simple lines repre- 
sent the soft-band emission (0.25-1.5 keV) while 
the lines with symbols represent the hard X-ray 
emission (1.5-7 keV). 
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Fig. 19.— Variation of soft-band (0.25-1.5 keV) 
X-ray luminosity with time in the WR phase for 
three sets of 40 M Q 2D models, assuming Solar 
abundances. Solid line: MM2003 with no rota- 
tion, dashed line: MM2003 with rotation, dotted 
line: STARS. In each case the simple lines repre- 
sent the model without thermal conduction, while 
the lines with symbols represent the models with 
conduction. 



of the soft X-ray emission from all of the 40 M Q 
2D models, assuming Solar abundances: MM2003 
with and without stellar rotation, and STARS, 
both with and without thermal conduction. Un- 
fortunately, in all cases, the hot gas moves off the 
2D computational grid (radius 10 pc) after a few 
tens of thousands of years, and this is what causes 
the reduction in the soft X-ray luminosity. How- 
ever, we comment that this is a larger size scale 
than either of the X-ray observed WR wind bub- 
bles. Our results are very varied. The STARS 
model, in which, as commented above, the shell 
of RSG material rapidly breaks up into clumps, 
has the lowest soft X-ray luminosity by an or- 
der of magnitude. In the MM2003 models, al- 
though the swept-up shells suffer instabilities, they 
do not break up into clumps before they move 
off the grid. The differences between the mod- 
els with and without thermal conduction should 
be attributed mainly to differences in the position 
(and volume) of the wind interaction region due 
to the effect of thermal conduction on the pres- 
sure in the main-sequence hot bubble, and only 
partly due to the local effect of thermal conduc- 
tion on the interaction region between the hot WR 
wind and the swept-up RSG material (see Fig. HT]) . 
Only three of our models attain X-ray luminosities 
above 10 33 erg s _1 : the MM2003 model without 
stellar rotation for the case without conduction, 
and the MM2003 models with stellar rotation both 
with and without conduction. Recall that the esti- 
mated total 0.2-1.5 keV band luminosity of S308 
is 1.2 x 10 34 erg s _1 (see Table [2| from analysis o f 
the NW quadrant of this object (|Chu et alJl2003h . 
although a recent study of the full set of observa- 
tions of S 308 suggests that the total luminosity is 
actually ~ 5 x 10 33 erg s -1 (see our companion 
paper Toala et al., in preparation). 

6. Discussion 

6.1. Comments on the stellar evolution 
models and their consequences for the 
circumstellar medium 

The stellar evolution tracks shown in Figure |4] 
indicate that none of the MM2003 models, 40 or 
60 Mq, with or without stellar rotation actually 
become red supergiants, since the minimum stel- 
lar surface temperature attained by these models 
is above 4800 K. Rather, they become yellow su- 
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Table 3 

Mass loss and duration of yellow supergiant and red supergiant stages 



Model 


M FMS a 


-WlWR b 


t YR c 


AM YR d 


tYB° 


AM YB f 


£rg s 


AM RG h 




M 


Mq 


yr 




yr 




yr 


Mq 


MM2003 


















40 M N' 


36.5 


15.7 


12100 


0.765 


391500 


18.40 






40 M R 


32.8 


22.2 


3420 


0.335 


54200 


8.00 






60 M N 


51.5 


26.7 


4410 


1.050 


10300 


3.21 






60 M Q R 


44.4 


32.2 














STARS 


















40 M N 


35.5 


16.7 


785 


0.042 


3210 


0.282 


77600 


15.8 


45 M Q N 


39.2 


18.6 


724 


0.069 


2690 


0.455 


31100 


16.1 


50 M Q N 


42.9 


20.4 


718 


0.119 


3400 


0.896 


16400 


16.1 


55 M Cl) N 


46.4 


22.5 


802 


0.244 


6570 


3.310 


8320 


13.3 


60 M Q N 


48.9 


23.8 


809 


0.325 


7730 


4.210 


6020 


12.0 



a Mass at the end of the main-sequence stage. 
b Mass at the beginning of the Wolf-Rayet stage. 
c Time spent as a yellow supergiant heading redwards. 
d Mass lost as a yellow supergiant heading redwards. 
c Timc spent as a yellow supergiant heading bluewards. 
f Mass lost as a yellow supergiant heading bluewards. 
s Time spent as a red giant. 
h Mass lost as a red giant. 

'N indicates models without stellar rotation, R indicates models with an initial rotation 
rate of 300 km s _1 at the stellar equator. 
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pergiants. The STARS models depicted in Fig- 
ure [4] all become red supergiants. In Table |3] we 
list the duration and mass lost in the yellow and 
red supergiant phases of evolution, separating the 
yellow supergiants into those heading to the red 
and those heading back to the blue. We use the 
temperature range 480 < T P ff < 7500 K t o define 
the yellow supergiants ( Drout et al. 20091 ). 

The yellow supergiant stage is a particularly in- 
teresting stage of stellar evolution, since statistical 
surveys of nea rby galaxies sugge st that it should 
be short lived (|Drout et al.ll2009h . For the STARS 
stellar evolution models, this is indeed the case, 
with the yellow supergiant phase lasting at most 
a few thousand years. Furthermore, the most im- 
portant mass loss occurs during the red supergiant 
phase, when the stellar wind velocities are lowest 
(< 30 km s _1 , although mass loss in the hotter, 
bluer phase becomes more important as the initial 
stellar mass increases). For the MM2003 models, 
on the other hand, there is no red supergiant stage, 
and the yellow supergiant stage lasts tens, or even 
hundreds, of thousands of years. For the MM2003 
40 Mq model without rotation, the most impor- 
tant mass loss occurs as a yellow supergiant head- 
ing bluewards. However, for the other MM2003 
models, the most important mass loss does not 
even occur in the yellow supergiant phase, but in 
a hotter, bluer evolutionary state. 

The different periods of mass loss have conse- 
quences for the circumstellar medium around the 
massive star. Intense mass loss in a short-lived 
red supergiant phase will result in a slow-moving 
dense shell of material close to the star; mass-loss 
spread out over hundreds of thousands of years 
in a yellow supergiant phase will lead to an ex- 
tended region of expanding circumstellar medium. 
Mass lost predominantly in hotter, bluer evolu- 
tionary stages may form fast-moving shells of ma- 
terial. When the fast WR wind interacts with 
the circumstellar medium, different structures will 
be produced depending on where the bulk of the 
circumstellar mass is and h ow fast it is moving 
( Garria-Segura et al. 1996al lb). 

From our results, we see that individual dense 
clumps are formed when the fast wind interacts 
with dense, slow-moving material close to the star, 
such as is the case for the STARS 40 and 60 M Q 
models (see Figs. [8] and [12]). This is because the 
linear thin-shell instabilities have time to break 



up the shell before it has traveled too far from 
the star. The clump densities and temperatures 
are such that the clumps will be photoionized and 
we would expect Ha emission in this case. When 
the circumstellar medium is moving more quickly, 
then the main fast wind-slow wind interaction oc- 
curs further from the star and, due to the lower 
densities (because of the geometrical dilution), the 
instabilities take longer to break up the swept-up 
shell. This is what we see in our MM2003 mod- 
els (Figs. |H [TUl [T3j) . In this case, the densities of 
the clumps and filaments are lower and their tem- 
peratures are higher, so we would expect neutral 
material only at the beginning of the interaction, 
when the shell is densest. 

6.2. Direct comparision with individual 
nebulae 

In Table [2] we present a summary of the ob- 
served and derived properties of three WR bubbles 
S308, NGC6888, and RCW58 and their central 
stars WR6, WR136, and WR40. In this section, 
we compare these properties and other details with 
the results of our numerical simulations. 

6.2.1. S308 

The wind-blown bubble S 308 around the WN4 
star WR6 is surprisingly spherical in [O III] im- 
ages, apart from a "blowout" region to the north- 
west. The Ha emission fr om this object shows 
a filamentary structure, and Gruendl et al. ( 2000f ) 
classified it as a Type II bubble. The optical ra- 
dius of S 308 is 9± 2 pc at a distance of 1.5±0.3 kpc 
(|Chu et al] 120031) . Our figures in § 14.21 were cho- 



sen such that the main features of the swept-up 
shell are roughly at this distance. It is important 
to point out that we are not trying to fit the wind 
parameters from the central star to the observed 
values, we are only choosing hydrodynamic results 
that approximately reproduce the morphology us- 
ing the wind parameters obtained from the most 
appropriate of our stellar models. 

The stellar evolution tracks shown in Figure 2] 
suggest that the central star, WR6, would be well 
fit by any one of the MM2003 40 M models with 
or with out stellar rotation, o r by a STARS 45 M 



model (Ham ann et al. I l2006h . Our numerical sim- 
ulations show that the STARS 45 M Q model (not 
shown, but very similar morphologically to the 
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0.6 0.8 
E [keV] 

Fig. 20.— X-ray spectra for the 40 M© 2D 
MM2003 model without rotation but with thermal 
conduction, corresponding to the DEM of Fig. [T7] 
Solid line: Modified, nitrogen-rich abundances, 
dashed line: Solar abundances, dotted line: Sin- 
gle temperature model (T = 3.6 x 10 6 K) with 
nitrogen-rich abundances. All spectra take into 
account absorption by a neutral column density 
of 2Vh = 1.1 x 10 21 cm~ 2 and are convolved with 
the EPIC/pn instrumental response matrices. 



40 M© model) would produce a Type III bubble 
(Ha clumps well separated from an [O III] shell). 
On the other hand, the MM2003 40 M© model 
without rotation would appear to produce struc- 
tures where dense, photoionized gas exists just in- 
side a warmer, collisionally excited (i.e., shocked) 
region, which could be interpreted as a Type II 
bubble. Model MM2003 with rotation breaks up 
the shell too far from the star and at the point 
depicted in Figure I10l (when the wind- wind inter- 
action is at 8 pc), the Ha emission and [O III] 
emission would be coincident. 

Spectroscopic analysis of S 308 shows the 
nebula to haye enh anced nitrogen abundances 
(jEsteban et al.lll992l) . Our simple estimation of 



the evolution of the nebular abundances with time 
(see Fig. H6| . made by time-averaging the total 
ejected mass of the different elements in the post- 
main-sequence phase, where the relative abun- 
dances are given by those at the stellar surface 
at a given time, suggests that both the STARS 
model and the MM2003 40 M© model with ro- 
tation could reproduce the observed abundances. 
However, the MM2003 model without rotation 
would not achieve the nebular abundances until 
the very end of its evolution. 

Another test of the models is to compare the 
predicted X-ray luminosity and spectrum with ob- 
servations. As can be seen from Figure \17\ our 
simulations produce large amounts of gas with 
temperatures T < 10 6 K, whereas observations are 
generally fit by a single temperature component, 
which is used to esti mate the to t al una bsorbed lu- 
minosity. For S 308. Ichu et al. ( 2003 ) estimate a 
total X-ray luminosity in the 0.25-1.5 keV band of 
Lx = 1.2± 0.5 x 10 34 erg s _1 by extrapolating their 
single temperature component (T ~ 1.1 X 10 6 K) 
results from the northwest quadra nt assuming a 
spherical bubble. IChu et al. (|2003h find that the 
spectrum is well fit using enhanced nitrogen ab un- 
dances, as determined bv lEsteban et al. 1 f)l992h for 
an absorbing column density of neutral hydrogen 
of 7V H = 1.1 x 10 21 cm" 2 (see Table [2]). From 
Figure US we see that the STARS 40 M© model 
completely fails to reproduce the observed X-ray 
luminosity, while the MM2003 model without ro- 
tation achieves this luminosity when thermal con- 
duction is not included, and the MM2003 model 
with rotation both with and without thermal con- 
duction produces sufficient X-ray luminosity in the 
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soft energy band. We remark that the results in 
Figure [19] were calculated for Solar abundances 
and that if the enhanced abundances are used, 
the luminosities are slightly different, depending 
on the dominant temperature. For example, the 
soft-band X-ray luminosity corresponding to the 
40 M Q MM2003 model with rotation and ther- 
mal conduction, whose differential emission mea- 
sure is shown in the lower panel of Figure HZ] 
is 8.2 x 10 33 erg s _1 for Solar abundances, but 
4.8 x 10 3 3 erg s -1 when we use the same abun- 



1.8 kpc ((van der Hucht et all 119881) . It is ellip- 



tical i n shape and was classified bv lGruendl et al 



dances as IChu et al.l (|20031 ). 

We also calculate the predicted X-ray spectra 
of our models. We note that none of our mod- 
els have been particularly tailored to S 308 (for 
example, we have not tried to match the stellar 
wind velocity to that of the central star WR6). 
The main difference between our models and the 
reported observations is that the simulations con- 
tain a range of temperatures that contribute more 
or less equally to the total spectrum, whereas 
Chu et al. [ (l2003h find that the observed spectrum 
can be well fit by a single temperature compo- 
nent at 1.1 x 10 6 K. In Figure [20] we show the 
simulated spectra for the 40 M Q MM2003 model 
without stellar rotation and with thermal conduc- 
tion, where we have considered Solar abundances 
and als o the modified, n itrogen-rich abundances 
used by IChu et al.1 (|2003f ). The assumed distance 
is 1.5 kpc and the absorption column density is 
iVn = 1-1 x 10 21 cm~ 2 , and we also convolved our 
model spectra with the instrumental response files 
appropriate for EPIC /pn on XMM-Newton. From 
this figure we see that the abundances make a 
large difference to the appearance of the spectrum. 
This is principally due to the reduced abundance 
of iron, which is responsible for the lines at around 
0.72 keV. The model spectra are very similar to a 
single temperature model with T ~ 3.6 x 10 6 K, 
which is hotter than that derived for S308. It is 
possible that non-equilibrium ionization effects are 
important in this nebula. The DEM for this simu- 
lation is shown in Figure [17] and we see that most 
of the gas actually has temperatures below 10 6 K 
but the emission from this gas is completely ab- 
sorbed by neutral hydrogen along the line of sight. 

6.2.2. NGC6888 

NGC 6888 is the wind-blown bubble around 
the WN6 star WR136, located at a distance of 



(2000) as having Type II morphology along the 
major axis but Type IV (no Ha emission) along 
the minor axis. The radius, determined from the 
[O III] image, is ~ 4.5 pc. Morphologically, this 
bubble is more like one of the MM2003 models 
than the STARS models, although the elliptical 
shape suggests that the mass loss is non-spherical. 

The central star of NGC 6888 sits well below 
any of our evolutionary tracks in the HR diagrams 
of Figure [4] However, the minimum mass for a 
star to become a Wolf-Rayet star is 37 Mq for 
MM2003 models without rota tion but 22 Mp, f or 



MM2 003 models with rotation ([Mevnet & Maeder 
2003). The nitrogen-to-oxygen abundance for this 
nebula shows a high level of nitrogen enrichment, 
which is also consistent with the stellar evolution 
models with rotation. 

NGC 6888 was detected i n X-rays by R OSAT 
(IWrigge et al.lll994l). ASCA (IWrigge et al.ll2005l) . 



and SUZAKU ([zhekov fc Parkll201ll) . Like S308, 
it has a limb-brightened X-ray surface brightness 
profile and a soft X-ray spectrum, which can be fit 
with two-temperature emission models where the 
principal component has T = 1.3 x 10 6 K and the 
second component has T = 8.8 x 10 6 K for an ab- 
sorption column density jVh = 3.1 x 10 21 cm~ 2 . 
The total X-ray luminosity in the 0.4-2.4 keV 
band is ~ 3.0 x 10 3 4 erg s -1 , where Sol ar metallic- 
ities were assumed ( Wrigge et al.ll2005 ). Although 
none of our simulations is consistent with the stel- 
lar parameters of this object, we remark that our 
models produce a range of temperatures in the X- 
ray emitting gas (see Fig. ITT|) , which could easily 
be interpreted as two main components. 



6.2.3. RCW58 

RCW58 is the wind-blown bubble around the 
WN8 star WR40. The optical emission is com- 
posed primarily of radial Ha filaments and clumps 
and the much fainter [O III] emission is smo other 



and more extended. Gruendl et al.l (|2000() clas- 



sify this object as having Type III morphology. 
The Ha filaments form a roughly elliptical struc- 
ture of d imension 6 x 8 pc, assuming a distance 
of 3 kpc ||Crrulll982t iHerald et ai]|200lh . which is 
roughly half the diameter of S 308. The morphol- 
ogy is similar to that given by the STARS models 
in Figures 151 and [T2l when the RSG shell has been 
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totally disrrupted. 

The stellar evolution tracks shown in Figure |4] 
suggest that the central star would be fit by ei- 
ther of the 60 Mq MM2003 models but the 60 M Q 
STARS is underluminous by an order of magni- 
tude and a higher initial mass model is suggested. 
The nitrogen-to-oxygen abundance determined for 
RCW 58 indicates material expelled from the star 
early on in the post-main-sequence evolution, since 
the nitrogen enhancement is only moderate. This 
could be due to the measurements being predom- 
inantly of clump material. 

RCW 58 was observed i n X-rays by XMM- 



Newton but not detected ( Gosset et al. 2005) 



Our Figures [T7] and \T§\ suggest that this could 
be because morphologies such as those produced 
by our STARS simulations do not produce much 
gas around T ~ 10 6 K. This seems to be because 
most of the mass of the RSG stage is bound up in 
the clumps. Clump material is photoevaporated 
by the WR ionizing flux and this material inter- 
acts with the fast wind. The resulting amount of 
material at ~ 10 6 K is rather small. Another con- 
tributing factor to the non-detection of RCW 58 is 
the fact that it is about twice as distant as S 308 
and hence the absorption column density, Nn, to 
this object is higher, which will have a large effect 
on the detection of soft X-ray emission. 

7. Summary and conclusions 

We have carried out numerical simulations in 
one and two dimensions of the evolution of the 
interstellar and circumstellar bubbles that form 
around massive stars during their lives. We find 
that the structures that form around the stars in 
the final stages of their lives are highly depen- 
dent on the details of the mass loss in the post- 
main-sequence stages and that this can be used 
to discriminate between different stellar evolution 
models. We have compared resu lts obtained with 
the p ublicly availa ble STARS (lEldridge fc Tout 
2004) and MM2003 dMevnet fc Maederll2003l) steb 



lar evolution models with observations of wind- 
blown bubbles around three Galactic Wolf-Rayet 
stars: S308, NGC6888 and RCW 58. In particu- 
lar we find that : 

1. We confirm the findings of previous authors 
that the interaction of a fast Wolf-Rayet 
wind with very slow moving RSG material 



results in the break up of the swept-up shell 
of material into numerous clumps and fila- 
ments, while the interaction of the fast wind 
with faster moving YSG or LBV material 
leads to corrugated shells of swept-up mate- 
rial several parsecs from the star. 

2. When thermal conduction is included in the 
simulations, the main differences occur in 
the main-sequence stage, since the timescale 
here is sufficient for the conduction pro- 
cess to occur. For the saturated conduction 
model that we use in this paper, conduction 
is most important at the edge of the bub- 
ble and the increased cooling rate that re- 
sults lowers the pressure in the hot bubble. 
Thus bubbles with thermal conduction ra- 
diate more energy and are slightly smaller 
than their counterparts without conduction. 
In the later stages of stellar evolution, con- 
duction does not affect the structures formed 
when the Wolf-Rayet wind interacts with the 
slow wind. 

3. We used stellar wind parameters from three 
sets of publicly available stellar evolution 
models for a range of initial stellar masses 
and find that the full range of morphologies 
are reproduced. However, other tests of the 
models are comparisons with optical spec- 
troscopy (e.g., the chemical abundances in 
the nebulae) and X-ray observations, both 
luminosities and spectra. 

4. We find that the STARS stellar evolution 
models, both 40 and 60 Mq, always result in 
the complete break up of the swept-up RSG 
material into clumps and filaments early on 
in Wolf-Rayet phase, since the slow wind ve- 
locities become so low (< 15 km s _1 ) that a 
large amount of dense material remains close 
to the star before the onset of the fast wind. 

5. None of the MM2003 evolution models, both 
with and without stellar rotation, result in 
clumpy wind-blown bubbles (for wind-wind 
interaction radius R ^ 10 pc) because the 
long duration of the YSG phase and faster 
wind velocity (~ 30 km s _1 ) result in a 
more extended distribution of dense mate- 
rial, which only becomes swept up into a 
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dense shell at several parsecs from the cen- 
tral star. 

The X-ray luminosities of bubbles with 
clumps and filaments are two orders of mag- 
nitude less than those with swept-up shells 
(10 32 erg s _1 as opposed to 10 34 erg s -1 ) and 
the gas is low temperature, hence the X-rays 
will be very soft. This could be the reason 
that RCW58 has not been detected in X- 



whereas S308 and NGC 
-l 



> 



rays 

10 34 erg s" 1 ) have been detected. Models 
with thermal conduction have slightly dif- 
ferent soft-band luminosities to their coun- 
terparts without conduction but there is no 
systematic trend. 

7. Observationally determined nitrogen-to- 
oxygen ratios reveal that NGC 6888 and 
S 308 have very enhanced nitrogen abun- 
dances compared to Solar but RCW58 is 
only moderately enriched. Of our stel- 
lar evolution models, the STARS models 
give enhanced abundances from the onset of 
the Wolf-Rayet phase, the MM2003 models 
without rotation give enhanced abundances 
towards the end of the Wolf-Rayet phase, 
and the MM2003 models with rotation give 
enhanced abundances even during the main- 
sequence stage due to mixing of processed 
material from the stellar interior to the sur- 
face. The models with rotation are unable to 
explain moderate abundances such as those 
determined for RCW 58. 

8. The observed X-ray spectra of S308 and 
NGC 6888 are very soft and as such must 
be strongly affected by interstellar absorp- 
tion. These spectra can be well fit by one, 
or at most two, components, with the main 
component having temperature T ~ 1.1— 
1.3 x 10 6 K. Our simulations, on the other 
hand, produce gas with a wide range of tem- 
peratures, which all contribute to the ob- 
served spectrum. 

Finally, we remark that none of the stellar evo- 
lution models we have tested can simultaneously 
reproduce all of the features of Wolf-Rayet wind- 
blown bubbles, that is morphology, abundances, 
X-ray luminosity and spectral characteristics. 



SJA and JAT acknowledge financial support 
from DGAPA, UNAM through project PAPIIT 
IN100309. JAT thanks CONACyT, Mexico for a 
scholarship. We thank Will Henney for a critical 
reading of the manuscript. We would like to thank 
the anonymous referee for constructive comments, 
which improved the presentation of this paper. 



37 



A. Details of the numerical method 



In this appendix we include more details of the numerical method for the interested reader 
A.l. Equations 

In spherical symmetry, the hydrodynamic equations in conservation form are 

dp 1 dr 2 pu 



dt r 2 dr 
dpu 1 dr 2 (p + pu 2 ) 2p 



dt r 2 dr 
de 1 dr 2 u(e + p) 



<7m, (Al) 

(A2) 



G - L + q e , (A3) 



dt r 2 dr 

where p, u, p are the mass density, radial velocity and pressure respectively, e is the total energy, defined by 

P 1 _„ 2 



- 7 -l ' 2^ (A4) 

G and L are the heating and cooling rates, and 7 is the ratio of specific heats. The stellar wind is injected 
as a thermal wind, with mass-loss rate q m — 3M w /47ri?^,, and thermal energy q c — 3i w /47ri?^, where M w is 
the stellar wind mass-loss rate, L w — M w V w 2 /2 is the stellar wind mech anical luminosity, V w is the stellar 
wind velocity, and R w is the radius of the stellar wind injection volume ( Chevalier fc Clegel [r985). Outside 
of the stellar wind injection region, i.e, r > i? w , we have q m = q c = 0. 

In our spherically symmetric simulations, we choose R w = 30dr, where dr is the grid cell size. In all our 
ID simulations, the cell size is set to be dr = 0.002 pc. The stellar wind injection region is small enough 
such that the stellar wind reaches its terminal velocity well before it shocks. This method of injecting the 
stellar wind is extremely robust to variations in the stellar wind parameters and does not result in spurious 
shock waves propagating back towards the center of symmetry. 

In addition to the usual hydrodynamic equations, we use advection equations for the neutral and ionized 
density, which couple the hydrodynamics to the radiative transfer, for example 

dpnx , 1 d / 2 



dt 



- 1 — (r 2 pun x )=0, (A5) 



where nx is the ionized or neutral number density. The ionized hydrogen fraction is updated using the 
current photoionization, collisional ionization and recombination rates in the computational cell. 

dy\ 

= Hed-^n + Pi-xVn - ™e-Ri?/i , (A6) 

where y n and yi are the neutral and ionized hydrogen fractions, n e is the electron density, Ci_i is the collisional 
ionization rate of neutral hydrogen and i?i is the radiative recombination rate of ionized hydrogen. Both 
Ci_i and i?i are functions of temperature only and we use analytic fits for both rates. The photoionization 
rate in each cell, Pj_i, is obtained from the radiative transfer procedure. 
In cylindrical symmetry, the conservation equations in the r-z plane are 

dp 1 drpu r dpu z 

9m, (A7) 



dt r dr dz 

dpu r 1 dr(p + pu 2 ) dpu r u z p 

dt r dr dz r ' 

dpu z 1 drpu r u z d(p + pup _ Q 

dt r dr dz 



(A8) 
(A9) 



38 



de + 1 dm, (e + p) + &x.(e + p) = g _ £ + 

at r or dz 
where « r , w z are the radial and azimuthal velocities, respectively, p is the pressure, and e is the total energy, 
defined by 

e =-P- + ip(«? + u2). (All) 
7 — 1 z 

The r and z directions are dealt with using operator splitting. 

A. 2. Heating and cooling 

The energy equation includes radiative cooling and heating due to the absorption of stellar radiation. 
In photoionized gas, the radiative cooling term should take into account the fact that hydrogen is fully 
ionized even when the gas temperature is only 10 4 K, and so collisional ionization of hydrogen is not an 
important cooling pro cess in an H II regi on (see Fig. IzTI) . We generate a cooling curve using the Cloudy 



photoionization code ([Ferland et al.l 119981 ) tailored to the ionizing source, which we tabulate and use as a 
look-up table during the calculation. Cloudy uses the most up-to-date and complete set of atomic data in 
the astrophysical literature and is constantly updated. 

The heating term at a given point in the numerical grid depends on the photoionization rate at that 
point and the stellar effective temperature, where the photoionization rate is calculated using the radiative 
transfer procedure. 

The heating and cooling are source terms in the energy equation, which is solved together with the 
hydrodynamic equations. If thermal conduction is included in the simulations, this is treated as an update 
to the energy equation using operator splitting after the hydrodynamic equations have been solved. 

A. 3. Radiative transfer 



In this work, we use the method of short characteristics, described in detail by iRaga et al. (|l999h . This 



method consists of finding the column density from a source point to any point on the grid by first calculating 
the column density at all intermediate points, beginning with those closest to the source. The column density 
at point P is therefore given in terms of the column density at a neighbouring point C by the equation 

N P =N C + ln ofi , (A12) 

where / is the path length from C to P and no,c is the neutral hydrogen density at point C . This method can 
be used both for the direct radiation from point sources and also for the diffuse radiation from a distributed 
source. In the spherically symmetric case, for direct radiation from a point source, the procedure is trivial. 
For cylindrical symmetry, the calculation of the path lengths I is slightly more involved, since it is necessary 
to find the points where the characteristics cross the cell boundaries by linear interpolation. 

The photoionization rates are needed to calculate the ionization state and heating rates of the gas. The 
column densities are therefore calculated at the beginning of every timestep and these are used to find 
the optical depths to each grid cell center and hence the photoionization rates at those positions. The 
update to the ionized (and neutral) hydrogen fractions is done by operating splitting after the solution of 
the hydrodynamic equations. 
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